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SADDLEPOINT  APPROXIMATION  AND  FIRST-ORDER  CORRECTION 


TERM  TO  THE  JOINT  PROBABILITY  DENSITY  FUNCTION 
OF  M  QUADRATIC  AND  LINEAR  FORMS  IN  K  GAUSSIAN 
RANDOM  VARIABLES  WITH  ARBITRARY  MEANS  AND  COVARIANCES 

INTRODUCTION 

When  K  normalized  Gaussian  random  variables  (RVs)  { g ( k ) }  are 
squared  and  summed,  the  resultant  z  is  called  a  chi-squared 
variate  with  K  degrees  of  freedom,  and  the  probability  density 
function  (PDF)  of  RV  z  is  available  in  a  closed  form  involving  an 
exponential.  If  constants  { c ( k ) }  are  added  to  each  of  the  RVs 
(g(k)}  before  squaring  and  summation,  the  PDF  of  the  resultant  z 
is  called  a  noncentral  chi-squared  variate  with  K  degrees  of 
freedom,  and  is  again  available  in  a  closed  form,  this  time 
involving  a  Bessel  function  and  an  exponential.  However, 
virtually  any  additional  complexity  beyond  this  case  results  in  a 
RV  z  for  which  the  corresponding  PDF  is  analytically  intractable. 

However,  in  these  one-dimensional  cases  of  RV  z,  the  moment 
generating  function  (MGF)  of  z  is  frequently  available  in  closed 
form,  and  a  numerical  technique  involving  fast  Fourier  transforms 
( FFTs )  can  be  efficiently  employed  to  get  numerous  quick  and 
accurate  values  for  the  PDF,  as  well  as  the  exceedance 
distribution  function,  at  arbitrary  points  of  interest,  whether 
near  the  mean  of  RV  z  or  on  the  tails  of  the  distribution  of  z 
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(references  1  through  5).  Thus,  the  one-dimensional  statistical 
problem  involving  quadratic  forms  of  Gaussian  RVs  is  in  good 
shape  numerically. 

The  situation  in  M  dimensions  is  much  more  difficult.  Even 
if  the  joint  M-dimensional  MGF  of  a  random  vector  (RV),  denoted 
by  column  vector  z  =  [ z ( 1 )  ...  z(M)]',  is  available  in  closed 
form,  its  inverse  M-dimensional  Laplace  or  Fourier  transform  back 
into  the  PDF  domain  cannot  be  accomplished  analytically,  except 
in  the  simplest  of  cases.  Also,  numerical  evaluation  of  the 
pertinent  M-dimensional  integral  for  the  joint  PDF  cannot  be  done 
accurately  for  M  greater  than  four  or  so.  These  conditions  force 
acceptance  of  an  approximation  to  the  M-dimensional  PDF  of  RV  z; 
also,  they  force  the  effort  to  be  concentrated  on  the  evaluation 
of  the  joint  PDF  at  very  few  points  in  M-dimensional  PDF  space, 
due  to  the  extensive  numerical  effort  and  execution  time  involved 
in  the  accurate  evaluation  of  multiple  integrals. 

The  M-dimensional  PDF  approximation  adopted  here  is  that 
obtained  via  the  saddlepoint  (SP)  method,  with  a  first-order 
correction  term  (reference  6,  page  180).  The  saddlepoint 
approximation  (SPA)  is  accurate  on  the  tails  of  the  joint  PDF,  as 
well  as  near  the  mean  of  the  distribution.  For  its  evaluation, 
the  SPA  requires  the  calculation  of  some  partial  derivatives  of 
the  joint  MGF  up  through  fourth-order;  evaluation  of  these 
derivatives  will  consume  much  of  the  effort  in  this  report. 
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EXAMPLE  PROBLEMS 


CORRELATION  ESTIMATES 

Let  w  =  [w(l)  ...  w( K ) ] '  be  a  Kxl  real  Gaussian  RV  with  Kxl 

mean  vector  r  and  KxK  positive-definite  covariance  matrix  R;  that 
is, 

E { w}  =  r  ,  E{(w  -  r ) ( w  -  r)'}  =  R  ,  (1) 

and  E {  }  denotes  an  expectation.  An  autocorrelation  estimate  of 
sequence  w  at  delay  n  is  available  according  to 

K 

a(n)  =  )  |  w(k)  w(k  -  n)  for  n=0 : K-l  .  (2) 

k=n+l 

Suppose,  for  example,  that  only  the  correlation  estimates  at 
delays  n  =  0,  1,  3,  and  7  are  of  interest;  that  is,  M  =  4  and 
RV  z  has  components 

z(l)  =  a( 0 ) ,  z ( 2 )  =  a( 1 ) ,  z(3)  =  a(3),  z(4)  =  a(7)  .  (3) 

The  problem  of  interest  is  to  obtain  the  joint  PDF  of  RV  z  for 
arbitrary  sample  size  K  and  statistics  r  and  R. 

The  quantities  in  equations  (2)  and  (3)  can  be  written  as 
quadratic  forms 

z(m)  =  w'  P(m)  w  for  m=l:M  ,  (4) 

where,  for  example,  KxK  matrices 


3 


■10  0...  ■ 

■0100...  • 

0  10 

1 

10  10 

P(l)  = 

0  0  1 

• 

• 

• 

P(2)  =  § 

0  10  1 

• 

• 

• 

Matrix  P(2)  is  nonzero  only  on  the  super-  and  sub-diagonals 
numbered  1;  matrix  P(3)  is  nonzero  only  on  super-  and  sub¬ 
diagonals  3;  and  matrix  P(4)  is  nonzero  only  on  super-  and 
sub-diagonals  7. 


If  the  sample  mean  is  subtracted  from  data  sequence  w  prior 
to  calculation  of  correlation  estimates  (2),  the  quadratic  forms 
for  RV  z  in  equation  (4)  still  hold,  but  the  elements  of  the 
matrices  {P(m)}  for  m=l:M  in  equation  (5)  are  changed.  For 
example,  the  j,k  element  of  matrix  P(l)  is  now  8^  ”  1/K  instead 
of  where  8^  is  the  Kronecker  delta;  the  remaining  matrices 

{P(m)}  are  more  complicated,  but  each  element  in  matrix  P(m)  can 
be  evaluated  by  means  of  a  single  sum. 


If  the  correlation  estimates  are  to  be  unbiased,  additional 
scale  factors  are  required  in  equations  (2)  or  (3).  Again,  the 
quadratic  forms  in  equation  (4)  are  appropriate,  but  the  elements 
of  matrices  {P(m)}  require  additional  calculations  involving  the 
particular  scale  factors  adopted. 

Equation  (2)  involves  an  aperiodic  correlation  of  data  w. 

The  extension  to  cyclic  correlation  estimates  {c(n)}  can  also  be 
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formulated  in  terms  of  quadratic  forms  (4).  Consider  the  cyclic 
correlation  estimate  at  delay  n  =  1: 

z(2)  =  c ( 1 )  =  a( 1 )  +  w( 1 )  w(K)  ,  (6) 

where  the  added  term  represents  wraparound.  This  RV  immediately 
fits  equation  (4)  if  the  two  corner  elements  (upper  right  and 
lower  left)  in  matrix  P(2)  (equation  (5))  are  changed  from  0  to 
1.  Instead,  for  delay  n  =  2,  c(2)  requires  four  changes  in  its  P 
matrix;  namely,  from  0  to  1  of  the  four  elements  immediately 
bordering  the  two  corner  elements.  This  procedure  extends  to  any 
delay  n;  the  corresponding  P  matrix  for  cyclic  correlation 
estimate  c(n)  will  have  K  Is  on  the  n-th  super-  and  sub-diagonals 
and  their  wraparound  extensions. 

Cross-correlation  estimates  from  two  different-length  data 
sequences  u  and  v  can  be  written  in  the  form 

z(m)  =  u'  A(m)  v  for  m=l:M  ,  (7) 

where  RV  u  is  Jxl,  RV  v  is  Kxl,  and  matrix  A(m)  is  JxK.  By 
defining  augmented  (J+K)xl  RV  w  as  [ur  v' ] ' ,  equation  (7)  may  be 
reformulated  as 

z(m)  =  w'  P(m)  w  for  m=l:M  ,  (8) 

where  matrix  P(m)  is  (J+K)x(J+K)  for  m=l:M.  Thus,  cross¬ 
correlation  estimates  obtained  from  two  different  sequences  can 
also  be  expressed  as  quadratic  forms  of  a  concatenated  sequence. 
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FILTERED  SQUARED  DATA 


Suppose  data  w  in  equation  (1)  is  processed  as  follows: 

2 

x  =  A  w  ,  y(n)  =  x(n)  for  n=l:N  ,  z  =  B  y  ,  (9) 

where  matrix  A  is  NxK ,  matrix  B  is  MxN,  and  y  =  [ y ( 1 )  ...  y(N)]'. 
Thus,  Nxl  RV  x  is  a  filtered  (linearly  transformed)  version  of 
data  w,  which  is  then  squared  and  subjected  to  additional 
filtering,  resulting  in  Mxl  RV  z.  The  problem  is  to  determine 
the  joint  M-dimensional  PDF  of  RV  z. 

By  combining  the  operations  in  equation  (9),  the  component 
RVs  (z(m)}  of  z  can  be  expressed  as 

K 

z(m)  =  >  |  P(m;k,l)  w(k)  w(l)  for  m=l:M  ,  (10) 

k ,  1=1 

where  constants 

N 

P(m;k,l)  =  )  ]  B(m,n)  A(n,k)  A(n,l)  for  m=l:M;  k,l=l:K.  (11) 

n=l 

I' 

Thus,  the  RVs  {z(m)}  in  equations  (9)  and  (10)  can  again  be 
expressed  as  quadratic  forms  (4),  where  KxK  matrices 

P(m)  =  [P(m;k,l);  k,l=l:K]  for  m=l:M  (12) 

in  terms  of  its  elements  calculated  in  equation  (11).  That  is, 
the  classical  filter-square-filter  operation  is  basically  a 
problem  in  finding  the  joint  PDF  of  several  statistically 
dependent  quadratic  forms. 
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PROBLEM  FORMULATION 

QUADRATIC  AND  LINEAR  FORMS  OF  INTEREST 

The  formulations  above  all  resulted  in  purely  quadratic  forms 
for  RV  z  =  [z(l)  ...  z(M)]'.  More  generally,  interest  here  will 
be  concentrated  on  the  M  quadratic  and  linear  (QAL)  forms 

z(m)  =  w'  P(m)  w  +  p(m)'  w  +  q(m)  for  m=l:M  ,  (13) 

where  RV  w  is  Kxl ,  matrix  P(m)  is  KxK ,  vector  p(m)  is  Kxl,  and 
scalar  q(m)  is  lxl.  Also,  every  matrix  P(m)  for  m=l:M  is 
symmetric  without  loss  of  generality  (see  appendix  A). 

RV  w  is  presumed  to  have  joint  Gaussian  statistics  with  Kxl 
mean  vector  r  and  KxK  covariance  matrix  R,  as  in  equation  (1). 
Thus,  equation  (13)  exhibits  the  most  general  second-order  forms 
in  correlated  Gaussian  RVs  with  arbitrary  statistics.  Since  all 
M  components  of  RV  z  in  equation  (13)  utilize  the  same  Kxl  RV  w 
but  in  different  combinations,  these  M  components  {z(m)}  are 
statistically  dependent  on  each  other,  in  addition  to  being 
non-Gaussian .  These  complications  are  what  force  the  need  to 
resort  to  an  approximation  for  the  desired  joint  PDF  of  RV  z. 

There  are  five  types  of  input  information  required  to 
completely  specify  the  QAL  problem  posed  in  equation  (13).  They 
are:  the  Kxl  mean  vector  r  of  Kxl  RV  w,  the  KxK  covariance  matrix 
R  of  Kxl  RV  w,  the  M  matrices  {P(m)}  of  size  KxK,  the  M  vectors 
(p(m)}  of  size  Kxl,  and  the  M  scalars  (q(m)}  of  size  lxl. 
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CONDENSATION  OF  QUADRATIC  AND  LINEAR  PROBLEM 


The  Kxl  RV  w  can  be  expressed  in  terms  of  a  set  of  K 
normalized  RVs  g  =  [g(l)  ...  g(K)]',  which  have  zero  mean  and  an 
identity  covariance  matrix,  according  to 

w  =  r  +  S'  g  ,  E{g}  =  0  ,  E{g  g'}  =  I  ,  (14) 

where  R  =  S'  S.  For  example,  KxK  matrix  S  can  be  the  Cholesky 
decomposition  of  positive-definite  covariance  matrix  R.  Then,  by 
substitution  of  equation  (14)  into  equation  (13),  there  follows 

z(m)  =  g'  C(m)  g  +  v(m) '  g  +  c(m)  for  m=l:M  ,  (15) 

where 


C(m)  =  S  P(m)  S'  (symmetric  KxK ) " 

v(m)  =  S[p(m)  +  2  P(m)  r]  (Kxl) 

c(m)  =  q(m)  +  p(m)'  r  +  r'  P(m)  r  (lxl); 


for  m=l:M  . 


(16) 


Now,  RV  z  in  equation  (15)  depends  only  on  the  three  types  of 
fundamental  quantities  given  in  equation  (16).  This  condensation 
or  pre-processing  of  the  input  information  will  prove  very  useful 
later  when  the  desired  joint  statistics  ( M-dimensional  PDF )  of  RV 
z  are  derived. 


If  mean  vector  r  =  0  and  all  vectors  p(m)  =  0  for  m=l:M, 
then  all  vectors  v(m)  =  0  for  m=l:M,  and  equation  (15)  reduces  to 
z(m)  *  g'  C(m)  g  +  q(m)  for  m=l:M.  This  is  called  the  purely 
quadratic  case;  its  SPA  and  first-order  correction  term  to  the 
joint  PDF  of  z  is  much  simpler  than  for  the  general  QAL  problem. 
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MOMENT  GENERATING  FUNCTION  OF  QUADRATIC  AND  LINEAR  FORMS 
Let  Mxl  vector  X  have  components 

X  =  [ X( 1 )  ...  X(M)]'  .  (17) 

The  joint  MGF  of  M-dimensional  RV  z  in  equation  (15)  is  defined 
as 


fj(\)  =  E{exp(X'  z)}  =  E^exp 


M 

E 

Lm=l 


X(m)  z(m) 


=  E{ exp[ g'  D( X)  g  +  t ( X) '  g  +  u(X)]}  , 


where 


M 

D(X)  =  )  1  X(m)  C(m)  (symmetric  KxK)  , 

m=l 

M 

t(X)  =  YU  Mm)  v(m)  ( Kxl )  , 

m=l 

and 

M 

u(X)  =  YU  Mm)  c(m)  (lxl)  . 

m=l 


(18) 


(19) 

(20) 


(21) 


The  constant  quantities  (C(m)},  (v(m)},  and  (c(m)}  were  defined 
in  equation  (16)  for  m=l:M.  It  should  be  observed  that  matrix 
functions  D(X),  t(X),,and  u(X)  are  linear  in  the  components 
(X(m)}  of  Mxl  vector  X.  The  problem  of  interest  now  is  to 
evaluate  the  K-dimensional  statistical  average  in  equation  (18) 
in  order  to  determine  the  M-dimensional  MGF  p{\)  of  RV  z  in 
closed  form. 
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Recall  from  equation  (14)  that  Kxl  RVs  w  and  g  are  related  by 
a  linear  transformation.  Since  RV  w  was  presumed  to  have 
Gaussian  statistics,  RV  g  must  also  have  Gaussian  statistics.  In 
fact,  from  equation  (14),  Gaussian  RV  g  has  a  zero  mean  vector 
and  an  identity  covariance  matrix.  The  joint  PDF  of  RV  g  is  then 

p(g)  =  (2n)-K//^  exp(-  g'  g/2 )  for  all  g  ,  (22) 

where  g  =  [g(l)  ...  g(K)]f  is  a  K-dimensional  field  point.  The 

pertinent  K-fold  integral  representation  of  equation  (18)  is 

fj(\)  =  (2n)"K/2  J  dg  exp[-.5  g'  Q(X)  g  +  t(X)'  g  +  u(X))  = 


exp[|  t(X)'  Q  ( X)  1  t(X)  +  u(  X)  J 
[det(Q(X) )  ]H 

where  symmetric  KxK  matrix  Q(X)  is  defined  as 


(23) 


Q( X)  -  I  -  2  D ( X)  .  (24) 

The  K-fold  integral  in  equation  (23)  for  the  joint  MGF  //(X) 
converges  only  if  all  K  of  the  eigenvalues  of  matrix  Q(X)  are 
positive.  Equivalently,  all  K  of  the  eigenvalues  of  matrix  D(X), 
defined  in  equation  (19),  must  be  less  than  1/2.  This  eigenvalue 
restriction  establishes  a  boundary  on  allowed  values  of  vector  X 
in  the  M-dimensional  X  plane;  in  particular,  the  origin,  X  =  0, 
is  always  an  allowed  point,  that  is,  a  point  at  which  the  joint 
MGF  (23)  exists. 
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Although  equation  (23)  is  a  closed-form  expression  for  the 
joint  MGF  of  RV  z  in  equation  (15),  it  contains  numerous  branch 
points  and  overlapping  essential  singularities  in  the  complex  X 
plane,  which  make  it  impossible  to  obtain  the  corresponding  joint 
M-dimensional  PDF  analytically.  Furthermore,  for  M  large,  it  is 
not  possible  to  perform  a  numerical  M-dimensional  FFT.  Thus,  it 
is  necessary  to  resort  to  the  SPA  for  the  desired  PDF  in  M 
dimensions . 

The  corresponding  joint  cumulant  generating  function  (CGF)  to 
joint  MGF  (23)  is 

X(  X)  =  log  fj(\)  = 

-  -  |  log  det (Q( X) )  +  |  t(X)'  Q(X)_1  t(X)  +  u(X)  .  (25) 

In  order  to  utilize  the  SPA  and  its  first-order  correction  term, 
partial  derivatives  of  joint  CGF  X(X),  with  respect  to  its  M 
components  {X(m)}  up  through  the  fourth  order,  must  be 
determined.  One  important  feature  of  these  derivations  is  that 
KxK  symmetric  matrix 

M 

Q (  X )  =  1-2  D (  X)  =1-2  YZ,  Mm)  C(m)  (26) 

m=l 

is  linear  in  components  {X(m)}.  Here,  equations  (24)  and  (19) 
were  used. 


11/(12  blank) 


M— DIMENSIONAL  SADDLEPOINT  APPROXIMATION 


M— DIMENSIONAL  SADDLEPOINT  IN  X  DOMAIN 

Suppose  that  the  joint  PDF  of  Mxl  RV  z  defined  in  equation 
(15)  is  desired  to  be  evaluated  at  a  general  Mxl  field  point 
z  =  (z(l)  ...  z  ( M )  ] '  .  This  PDF  value  is  given  in  terms  of  the 
joint  MGF  fj(\)  by  the  M-dimensional  integral 

p(z)  -  — - — rr  [•••[  d\(  1 )  ...  dX(M)  exp[ -  X'  z]  /j(\)  = 

( i2n)n  J 

^1  lM 

=  — — jj  [**•[  dX(l)  ...  dX(M)  exp[X(X)  -  X'  z]  ,  (27) 
(i2n)n  J  J 

U1  UM 

where  contour  Cm  in  the  complex  X(m)  plane  goes  from  -i®  to  +i® 
and  stays  within  the  analytic  boundary  of  the  joint  MGF  //(X). 

The  SPA  consists  of  locating  these  contours  so  that  they  pass 
through  the  M-dimensional  SP  of  the  integrand  of  equation  (27), 
and  then  approximating  the  integrand  values  on  the  contours  by  a 
Gaussian  M-dimensional  mountain  in  the  neighborhood  of  the  peak 
at  the  SP.  Finally,  this  Gaussian  approximation  is  extended  to 
all  values  on  the  contours  of  integration,  for  which  the  modified 
M-dimensional  integral  is  capable  of  evaluation  in  closed  form. 

In  order  to  determine  the  SP  in  the  X  plane,  it  is  necessary 
to  find  the  location  of  the  minimum  of  the  real  quantity 

M 

X(  X)  -  X'  z  =  X( X)  -  X(m)  z (m)  (28) 

m=l 
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in  equation  (27)  for  a  real  X  vector.  Alternatively,  this  SP 
location  is  found  by  solving  the  M  simultaneous  nonlinear  real 
equations 


3X(_X) 
9X(  m ) 


z(m)  for  m=l:M 


{X(m)}  real  . 


(29) 


The  real  solution  X  =  X(z)  =  [X(l)  ...  X(M)]'  is  a  function  of 
the  particular  M-dimensional  field  point  z  of  interest.  If  this 
field  point  is  changed,  the  M  nonlinear  equations  (29)  must  be 
re-solved  for  the  new  SP. 


M-DIMENSIONAL  SADDLEPOINT  APPROXIMATION  TO  PDF  OF  z 


When  the  M-fold  integration  procedure  above  is  carried  out, 
the  resulting  SPA  to  the  joint  PDF  is  (reference  6) 


p(  z) 


A  A 

exp[x(X)  -  X'  z]  s  /_, 
(  2 n ) M/2  (det(H(  X) )  ]**  0 


(30) 


where  H(X)  is  the  MxM  symmetric  Hessian  matrix  of  second-order 
partial  derivatives  of  the  joint  CGF : 


H(  X) 


92X( X) 
3X(m)  3X(m) 


m,m=l:M  . 


(31) 


The  function  Pq(z)  is  denoted  as  SPAO,  meaning  the  zeroth-order 
SPA  to  the  joint  PDF  p(z);  this  nomenclature  distinguishes  it 
from  some  further  approximations  to  the  joint  PDF  p(z)  that  will 
employ  a  first-order  correction  term. 
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It  is  useful  to  define  a  Gradient  vector  as 


G(X) 


raxiAI  ...  ax(\n  . 
Laxa)  3  x(  m  )  J  ' 


(32) 


then,  SP  equation  (29)  can  be  succinctly  expressed  as  G(X)  =  z. 

A  A 

If  the  solution  for  the  M-dimensional  SP  location  X  =  X(z)  is 
obtained  by  using  the  Newton-Raphson  search  procedure,  then  both 
the  Mxl  Gradient  vector  G(X)  and  the  MxM  Hessian  matrix  H(X)  will 
be  required  during  the  complete  search  procedure.  This 
necessitates  the  evaluation  of  first-  and  second-order  partial 
derivatives  of  the  joint  CGF,  as  indicated  in  equations  (31)  and 
(32). 


FIRST-ORDER  CORRECTION  TERM  TO  THE  SP  APPROXIMATION 


For  integers  j,k,l,m=l:M,  define  the  following  quantities, 

A  A 

which  are  evaluated  at  the  SP  X  =  X(z),  once  it  has  been 
determined  for  a  given  field  point  z: 


y  =  3X(X) 
Am  3X(m) 


32X(X) 


lm  3X(1)  3X(m) 


(33) 


3  X(X) 


vklm  "  3X(k)  3X(  1 )  3X(m) 


94X(X) 


jklm  “  3X( j )  3X( k )  3X(1)  3X(m) 


The  latter  two  quantities  do  not  need  to  be  evaluated  during  the 
search  for  the  SP,  but  only  need  to  be  evaluated  after  the  search 
has  been  completed.  Also,  define  the  two  symmetric  MxM  matrices 
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H  -  [XlB]  ,  T  -  a’1  -  tTlm]  . 
Finally,  define  the  three  constants 

°4  =  *  j5m  Xjklm  Tjk  Tlm  ' 


'3a 


r  2  Ex. 


8  Am  Am  *kl”  X“S  Tkl  V  *15  ' 


and 


'3b 


1 7  .2  I  Xi 


12  Am  Am  Aklm  *“2  V  Tli  V  ' 


where  the  sums  all  run  from  1  to  M. 


(34) 


(35) 


The  first-order  correction  to  the  SPAO  given  in  equation  (30) 
can  now  be  expressed  in  the  form  (reference  6,  page  180) 

P1(z)  s  pQ(z)  t1  +  ctl  »  (36) 

where  the  total  first-order  correction  term  is  defined  as 

ct  =  c4  +  c3a  +  c3b  '  (37) 

The  joint  PDF  approximation  p1(z)  in  equation  (36)  is  denoted  as 
SPAl,  meaning  the  first-order  SPA.  Its  computation  requires 
determination  of  the  SP  location,  as  well  as  third-  and 
fourth-order  information  about  the  partial  derivatives  of  the 

A  A 

joint  CGF  X(X)  at  the  SP  X  =  X(z). 
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MODIFIED  SADDLEPOINT  APPROXIMATIONS 


Consider  the  following  three  modified  SPAs: 
p1(z)  s  p0 ( z )  [1  +  c t ]  =  Pq ( z )  [1  +  cfc  +  0  cfc2  +  0  cfc3  + 

12  13 

pe(z)  s  Pq ( z )  exp(ct)  =  pQ ( z )  [1  +  cfc  +  j  cfc  +  g  cfc  + 
and 

1  +  ct//^  12  13 

pa(z)  s  P0(2)  1  -~c't72  =  P0(Z)  [1  +  ct  +  2  Ct  +  4  Ct  + 


]  , 

]  , 


]  . 
(38) 


Approximation  p^(z)  defined  in  equation  (36)  tacitly  employs  a 
zero  coefficient  for  the  second-order  correction  term  (SOCT)  ct2> 
This  coefficient  is  most  certainly  incorrect  because  there 
definitely  is  a  nonzero  SOCT;  however,  this  SOCT  is  not  known. 
Furthermore,  the  SOCT  c^2  would  require  knowledge  of  the  fifth- 
and  sixth-order  partial  derivatives  of  the  joint  CGF  X(M  at  the 
SP.  Since  there  are  M^  sixth-order  partial  derivatives,  a 
problem  arises  in  execution  time  and  storage  when  attempting  to 
calculate  these  latter  quantities. 


To  circumvent  the  lack  of  knowledge  and  computational 

limitations,  approximation  pg(z)  in  equation  (38)  has  been 

suggested  (reference  6,  page  180)  because  it  injects  the  SOCT 
2 

c^.  /2  instead  of  zero.  Again,  this  term  is  most  certainly 
incorrect;  however,  it  may  give  a  better  approximation  to  the 
true  joint  PDF  p(z)  than  either  of  the  approximations  Pq(z)  or 
P1(z) . 
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The  third  approximation,  p  (z),  in  equation  (38)  uses, 
instead,  a  rational  function  in  cfc,  which  has  the  same  power 
series  expansion  as  exp(c^)  through  second  order.  As  c^. 
increases,  approximation  pa(z)  becomes  greater  than  pg(z)  and 
would  tend  to  infinity  if  cfc  approaches  2;  however,  by  this  time, 
c^  could  no  longer  be  considered  a  correction  term,  but  in  fact, 
a  dominant  contributor. 

The  author  has  conducted  some  numerical  comparisons  of  the 
three  approximations  in  equation  (38)  for  some  cases  where  the 
exact  joint  PDF  p(z)  can  be  determined.  These  results  indicate 
that  both  Pe(z)  and  Pa(z)  generally  yield  worthwhile  improvements 
relative  to  p1(z),  which,  in  turn,  yields  worthwhile  improvements 
compared  to  Pq(z).  The  choice  between  pg(z)  and  pa(z)  varies 
from  example  to  example. 
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EVALUATION  OF  PARTIAL  DERIVATIVES  OF  JOINT  CGF  X(X) 


Evaluation  of  the  various  SPAs  depends  heavily  upon  the 
ability  to  obtain  the  partial  derivatives  of  the  joint  CGF  x(X) 
of  RV  z;  see  equations  (30)  through  (33).  This  joint  CGF  is 
repeated  from  equations  (25),  (26),  (20),  and  (21): 

X( X)  -  -  |  log  det ( Q( X) )  +  |  t(X)'  Q(X)-1  t(X)  +  u(X)  ,  (39) 

where 


M 

Q( X)  =1-2  D( X)  =  I  -  2  YZ  X(m)  C(m)  ( KxK )  ,  (40) 

m=l 


M 

t  ( X)  =  YU  X(m)  v<m)  (Kxl)  ,  (41) 

m=l 


M 

u( X)  =  YU  X(m)  c(m)  (lxl)  .  (42) 

m=l 


A  POSSIBLE  APPROACH  TO  THE  PARTIAL  DERIVATIVES 

From  equation  (25),  joint  CGF  X(X)  =  log  //(X);  therefore, 

3X(X) _ 1_  3//(X)  , 

3X(  m)  ”  /j(\)  3X(m)  * 

Also,  from  equation  (18), 

fj(\)  =  E{exp[X(l)  z(l)  +  •••  +  X(M)  z(M)  ] }  .  (44) 

There  follows  immediately 

=  E{  z(  m)  exp[  X(  1 )  z  ( 1 )  +  •••  X(M)  z(M)]}  .,  (45) 
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and 


32//(X) 
3X(1)  9X( m) 


E { z ( 1 )  z(m)  exp[ X( 1 )  z( 1 )  + 


X(M)  z  ( M )  ]  }  .  (46) 


Recall  from  defining  equation  (13)  and  condensed  version  (15) 
that  RV  z  is  quadratic  in  Gaussian  RVs;  therefore,  the  arguments 
of  the  exponentials  in  equations  (45)  and  (46)  are  quadratic  in 
Gaussian  RVs.  Similarly,  the  leading  multiplying  factor  z(m)  in 
equation  (45)  is  quadratic  in  Gaussian  RVs.  Therefore,  the 
multiple  integral  representing  expectation  (45)  can  certainly  be 
evaluated  in  closed  form,  in  a  similar  fashion,  the  leading 
factor  z ( 1 )  z(m)  in  equation  (46)  is  quartic  in  Gaussian  RVs, 
meaning  that  it  too  can  be  evaluated  in  closed  form.  The  same 
conclusion  holds  for  all  the  higher  order  partial  derivatives  of 
joint  CGF  X(X),  although  the  integral  evaluations  will  be 
considerably  more  tedious  to  carry  out. 


The  significance  of  this  observation  is  that  all  of  the 
required  information  for  obtaining  the  various  SPAs  is 
obtainable,  somehow,  in  closed  form.  The  best  route  for  getting 
this  information  may  not  be  by  means  of  expectations  (45)  and 
(46),  but  at  least  it  is  now  known  that  the  desired  information 
is  obtainable.  (An  example  of  this  route  to  the  partial  deriva¬ 
tives  of  the  joint  CGF  is  given  in  appendix  B;  some  interesting 
results  on  partial  derivatives  of  eigenvalues  are  also  provided.) 
However,  the  alternative  technique  described  below  is  much  more 
efficient  and  more  readily  supplies  the  required  higher  order 
joint  CGF  partial  derivatives  needed  for  the  various  SPAs. 
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FIRST-ORDER  PARTIAL  DERIVATIVES  OF  CGFl 


The  joint  CGF  x(X)  is  given  in  equation  (39).  It  is  composed 
of  three  additive  parts,  to  be  labeled  x^(X),  X2  ( M  ,  and  • 

The  partial  derivative  of  the  third  part,  Xg(X)  =  u(X),  with 
respect  to  X(m),  is  simply  c(m),  as  seen  from  equation  (42).  The 
higher  order  derivatives  of  Xj(X)  are  all  zero  because  {c(m)J  are 
constants  (see  equation  (16)). 

The  immediate  interest  in  this  subsection  is  in  the  first 
part,  CGFl,  of  the  complete  joint  CGF  X(X),  namely, 

XX(X)  =  -  j  log  det  Q(X)  ,  (47) 

where  the  symmetric  KxK  matrix  Q(X)  is  given  by  equation  (40). 

In  order  to  streamline  the  following  derivations,  a  number  of 
useful  matrix  properties  were  collected  in  appendix  C  and  will  be 
referred  to,  as  necessary. 

If  symmetric  KxK  matrix  D(X)  is  expanded  in  its  eigen- 
decomposition,  the  result  is 

M 

D ( X )  =  YU  X(m)  C(m)  =  V(X)  E(X)  v<x)'  '  V(X)  V( X) '  -  I  ,  (48) 
m=l 

where  KxK  matrix  V(X)  is  the  set  of  eigenvectors,  and  diagonal 
KxK  matrix  E(X)  is  the  set  of  eigenvalues  {e^(X)},  k=l:K.  There 
follows 

Q ( X)  =1-2  D( X)  =  V( X)  [1-2  E( X) ]  V ( X ) '  ,  (49) 
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because  matrix  V  has  a  unity  determinant.  Equation  (47)  yields 


K 

Xl(X)  -  -  |  log  det  Q(X)  =  -  4  log[l  -  2  ek(X)']  = 


k=l 


1  r — i  2P  /\\P  1  t — 1  2P  r-— 1  _  1  t — 1  2P  .  ,  ,,  ,p, 

1 S  R  5"  k<  ’  ■ 1 R  5~ £  e*u>p  -i£  r tt(D(X)  1 


=  4  trof]  ^  D(X)P1  =  \  tr  {-  log[  I  -  2  D(X)]}  , 


P 


lp=l 


(51) 


where  equations  (C-6)  and  (C-12)  were  used.  Now,  by  using 
equations  (C-18),  (48),  and  (49),  there  follows  the  desired 
first-order  partial  derivative 


3X1(X)  1 

3X(m)  “  2 


( 


tr  [I  -  2  D( X) ] 


-1 


-  3D( X) 1 
z  3X(m) J 


tr{Q(X)  1  C(m) }  (52) 


for  m=l:M.  The  symmetric  KxK  matrices  {C(m)}  are  given  in 
equation  ( 16 ) . 


SECOND-ORDER  PARTIAL  DERIVATIVES  OF  CGFl 

For  convenience  of  notation,  let  symmetric  KxK  matrix 
-1  M 

q(X)  =  Q(  X)  ,  Q(  X)  =  1-2  D(X)  =  YZ  X(m)  C(m)  .  (53) 

m=l 


Then,  by  using  equations  (C-15),  (48),  and  (49),  there  follows 


f^|  =  -  q(X)  |^j-  q  ( X)  =  2  q(X)  C(m)  q(X)  for  (54) 


This  enables  equation  (52)  to  be  written  in  the  compact  form 


3X1(X) 
3X(  m) 


tr { q ( X)  C(m) } 


tr{B(X,m)}  for  m=l:M  , 


(55) 


where  nonsymmetric  KxK  matrices  {B(X,m)},  m=l:M,  are  introduced 
for  future  use.  The  desired  second-order  partial  derivative  now 
follows  from  equations  (54)  and  (55)  as 


a2x1(X) 

9  X( 1 )  aX(m) 


=  2  tr{q(X)  C(l)  q(X)  C(m)}  =  2  tr{B(X,l)  B(Xfm)} 


for  l,m=l:M 


(56) 


THIRD-  AND  FOURTH-ORDER  PARTIAL  DERIVATIVES  OF  CGFl 

By  using  equations  (56)  and  (54),  it  immediately  follows 
that,  for  k,l,m=l:M, 


83X1(X) 


3 X (  k )  9X(  1 )  9X(m) 


8  tr {B ( X, k )  B ( X, 1 )  B( X,m) }  =  8  tr{Bk  B^  BmJ  , 

(57) 


where  the  shortcut  notation  B(X,m)  =  Bm  has  been  introduced. 

m 

With  this  notation,  the  final  quantity  of  interest  is 
34X1(X) 

3X( j )  3X(k)  3 X ( 1 )  9X(m)  =  16  tr{Bj  Bk  B1  Bm  + 

+  B_.  BR  Bm  Bj^  +  Bj  Bj^  Bk  Bm}  for  j,k,l,m=l:M  .  (58) 
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A  summary  of  the  notation  is  given  by 


Bm  =  B(X,m)  =  q( X)  C(m)  =  Q(X)  1  C(m)  =  [I  -  2  D(X)]  1  C(m)  (59) 

for  m=l:M,  with 
M 

D(X)  =  X(m)  C(m)  .  (60) 

m=l 

In  writing  expressions  (57)  and  (58),  advantage  has  been 
taken  of  symmetries  of  some  expressions  involved  in  matrices 
{Bm};  this  allowed  equation  (57)  to  be  condensed  into  a  single 
trace.  However,  the  three  traces  remaining  in  equation  (58)  are 
all  different  in  general,  and  no  further  reduction  is  possible  in 
the  number  of  terms  that  must  be  calculated. 


FIRST-  AND  SECOND-ORDER  PARTIAL  DERIVATIVES  OF  CGF2 


The  interest  is  now  centered  on  the  second  part,  CGF2 ,  of  the 
complete  joint  CGF  X(X)  in  equation  (39),  namely, 

j' 

1  M 
X2(X)  =  I  t  (  X)  '  q(t)  t  ( X)  ,  t(  X)  =  5m  Mm)  v(m)  .  (61) 

m=l 


The  pertinent  partial  derivatives  required  are  given  by  equation 


(54)  and 


3t(X) 

9X(m) 


v(m)  for  m=l:M  . 


(62) 


Application  to  equation  (61)  yields  the  first-order  partial 
derivative  of  CGF2  in  the  form 
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(63) 


8X2(  M 

9X(m) 


v(m)'  q  t  +  t'  B_a  t  for  m=l:M  , 

m  ^ 


using  an  obvious  shorthand  notation.  The  corresponding  second- 
order  partial  derivative  is  obtained  by  repeated  applications  of 
the  above  rules: 


32X2(X) 
9X( 1 )  9X( m) 

+  2  v(m) ' 


v(l) '  q  v(m)  +  2  v(l)'  Bm  q  t  + 

m 

qt+4t'B.B  qt  for  l,m=l:M  . 
1  m 


(64) 


It  should  be  noted  that  the  original  CGF2 ,  namely  quadratic 
form  X2(X)  in  equation  (61),  began  and  ended  with  Kxl  vector  t; 
this  is  called  a  (t-t)  type  of  term.  •  On  the  other  hand,  the 
first-order  partial  derivative  in  equation  (63)  involved  an 
additional  type  of  quadratic  form  starting  with  Kxl  vector  v(m) 
and  ending  with  vector  t;  this  is  called  a  (v-t)  type  of  term. 
Finally,  the  second-order  partial  derivative  in  equation  (64) 
involved  still  another  type  of  quadratic  form,  beginning  and 
ending  with  two  v  vectors;  this  is  called  a  (v-v)  type  of  term. 
At  this  point,  steady  state  is  reached;  that  is,  no  more 
additional  types  of  terms  are  generated  by  taking  additional 
higher  order  partial  derivatives  of  equation  (64).  However,  the 
numbers  of  each  type  of  term  do  increase,  and  the  complexity  of 
each  term  also  increases. 
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THIRD-ORDER  PARTIAL  DERIVATIVES  OF  CGF2 


Upon  taking  the  next  partial  derivative  of  equation  (64)  and 
combining  like  terms,  the  result  is 


33X2(X) 

3X(k)  3X(  1 )  3X( m)  = 

q  Ck  q  Cx  q  Cm  q  t  + 

k  m  1 
lkm 

(65) 


2  v/  q  C,  q  vm  +  4  v/  q  C,  q  C  q  t  +  8  t' 


m 

1 

m 

k 

1 

k 


for  k,l,m=l:M,  where  only  the  subscripts  have  been  indicated 
after  the  first  line.  There  are  31  =6  terms  of  type  (v-t),  but 
only  three  terms  of  types  (v-v)  and  (t-t).  The  reason  for  this 
is  that  the  transposes  of  half  of  the  (v-v)  and  (t-t)  (scalar) 
terms  can  be  shown  to  be  equal  to  those  displayed  in  equation 
(65);  these  common  values  have  been  combined  and  the  appropriate 
scale  factor  adjusted.  All  the  quadratic  forms  remaining  in 
equation  (65)  are  different  in  general;  no  further  reduction  in 
the  number  or  types  of  terms  is  possible. 
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FOURTH-ORDER  PARTIAL  DERIVATIVES  OF  CGF2 

The  fourth-order  partial  derivative  of  CGF2 ,  namely 
also  contains  only  the  (v-v),  (v-t),  and  (t-t)  terms.  In 
particular,  for  j,k,l,m=l:M,  the  (v-v)  terms  are 

4  v'.  q  C,  aC.qv  (66) 

j  ^  k  1  m 

and  12  permutations  of  its  subscripts.  The  (v-t)  terms  are 

e  q  ck  q  Cx  q  cm  ,  t  (67) 

and  24  permutations  of  its  subscripts.  The  (t-t)  terms  are 

8  t'  q  C.  q  Ck  q  Cj  q  Cm  q  t  (68) 

and  12  permutations  of  its  subscripts. 

There  are  41  =  24  terms  of  type  (v-t),  but  only  12  terms  of 
types  (v-v)  and  (t-t).  The  reason  is  identical  to  that  cited 
under  equation  (65),  namely,  the  equality  of  some  transposes  of 
scalar  quantities  involving  (v-v)  or  (t-t)  terms.  The  twelve 
quadratic  forms  remaining  in  equations  (66)  through  (68)  are 
different  in  general;  no  further  reduction  in  the  number  or  types 
of  terms  is  possible. 

It  should  be  noted  in  equations  (66)  through  (68)  that  a 
large  number  of  matrix  multiplications  are  involved,  especially 
in  equation  (68),  where  11  terms  are  involved.  However,  by 
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starting  at  one  end  of  equation  (68),  for  example,  all  the 
successive  multiplications  involve  a  vector  with  a  matrix,  which 
is  considerably  faster  than  for  full  KxK  matrices. 

Alternatively,  the  (Bm)  matrices  in  equation  (59)  can  be  computed 
once  and  stored  for  repeated  use  in  the  operations  above.  The 
danger  with  this  latter  approach  is  the  possibility  of  very  large 
storage  requirements,  especially  for  large  M  and/or  K. 

The  totality  of  partial  derivatives  required  to  compute  the 
first-order  correction  term  c^  to  the  SPA  in  equations  (33) 
through  (38)  is  given  in  equations  (55)  through  (58)  and 
equations  (63)  through  (68).  These  results  have  been  combined  in 
a  MATLAB  program  listed  in  appendix  D  and  entitled  quadlinspa, 
denoting  the  SPA  for  the  M  general  quadratic  and  linear  forms  of 
equation  (13).  In  the  special  case  of  purely  quadratic  forms 
(see  bottom  of  page  8),  an  alternative  MATLAB  program,  listed  in 
appendix  E  and  entitled  quadspa,  has  been  written.  Both  programs 
compute  the  three  SPAs  in  equation  (38)  as  well  as  the  standard 
SPAO  in  equation  (30). 
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SUMMARY 


The  saddlepoint  approximation  to  the  joint  M-dimensional 
probability  density  function  for  M  arbitrary  quadratic  and  linear 
forms  in  K  Gaussian  random  variables  with  arbitrary  means  and 
covariance  matrix  has  been  derived.  Also,  the  first-order 
correction  term  to  the  standard  saddlepoint  approximation  in  M 
dimensions  has  been  determined  and  used  to  form  several  different 
possible  approximations  to  the  desired  joint  probability  density 
function. 

The  determination  of  the  M-dimensional  saddlepoint  location, 
and  the  standard  saddlepoint  approximation  itself,  require 
evaluation  of  first-  and  second-order  partial  derivatives  of  the 
joint  cumulant  generating  function  at  arbitrary  points  in 
M-dimensional  space;  these  quantities  have  been  derived  in  closed 
form.  Also,  the  third-  and  fourth-order  partial  derivatives  of 
the  joint  cumulant  generating  function  have  been  derived  for  use 
in  calculating  the  first-order  correction  term  to  the  saddlepoint 
approximation  in  M  dimensions.  All  these  results  have  been 
combined  in  two  MATLAB  programs;  namely,  quadlinspa,  which 
handles  the  quadratic  and  linear  case,  and  quadspa,  which  handles 
the  purely  quadratic  case. 

Sometimes,  interest  is  centered  on  the  square  roots  of  the 
quadratic  and  linear  random  variables  {z(m)}  when  all  of  these 
quantities  are  nonnegative.  For  example,  in  some  signal 
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processing  applications,  the  square  roots  represent  envelope  or 
amplitude  quantities  of  interest,  while  the  {z(m)}  are  power 
quantities.  Letting  u(m)  =  z(m)'5  for  m=l:M,  the  joint 
probability  density  function  p2  of  the  random  variables  {u(m)}  at 
M-dimensional  field  point  u  =  [u^,...,u„]'  is  given  by 


PoU.,  ,  •  • 


'V 


(  2  21  -M 

-  plv-'vl  2 


U, 


u 


M 


(69) 


for  um  >  0  for  m=l:M.  Thus,  if  joint  probability  density  function 

P2  is  to  be  determined  at  arguments  u^,...^^  the  joint 

probability  density  function  p  of  random  vector  z  must  be 

2  2 

evaluated  at  arguments  u^,...,uM;  this  serves  as  the  field  point 

2  2 

z,  namely,  z  =  [u^, . . . ,uM] ' ,  for  the  procedures  detailed  above  in 
this  report.  More  generally,  this  procedure  can  be  extended  to 
nonlinear  transformations  u(m)  =  fm(z)  for  m=l:M,  provided  that 
the  right-hand  sides  {fm(z)}  do  not  generate  imaginary  numbers 
for  some  values  of  random  vector  z. 
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APPENDIX  A  -  A  PROPERTY  OF  QUADRATIC  FORMS 


Let  v  be  a  Kxl  vector  and  let  C  be  a  KxK  matrix,  not 
necessarily  symmetric.  The  symmetric  and  anti ( skew)-symmetric 
matrices  of  C  are  defined  as 

Cs  =  !<c  ♦  CM  ,  Ca  -  !<c  -  CM  . 

It  immediately  follows  that  C  =  C  +  C  and  C' 

5  a  a 

Now,  consider  the  quadratic  form 
f  =  v'  C  v  =  v'  (C  +  C  )  v  =  v'  C  v  for  any  v  .  (A-2) 

S3  S 

The  term  involving  C  is  zero,  as  can  be  seen  by  taking  its 
transpose  and  using  the  property  =  -C  ;  thus,  only  the 
symmetric  part  of  matrix  C  is  active  in  quadratic  form  f. 
Therefore,  when  a  quadratic  form  such  as  v'  C  v  is  encountered, 
the  matrix  C  may  be  presumed  symmetric  without  loss  of 
generality. 


( A— 1 ) 


Ca* 


A— l/( A— 2  blank) 


APPENDIX  B  -  DIRECT  EVALUATION  OF  EXPECTATIONS 


FIRST-ORDER  PARTIAL  DERIVATIVES  OF  CGF 

The  first-order  (FO)  partial  derivative  (PD)  of  joint  MGF 
fj{\)  is  given  by  the  expectation  in  equation  (45).  Also,  RV  z(m) 
is  given  in  equations  (14)  and  (15)  as 


z(m)  =  g'  C(m)  g  for  ra=l:M  , 


( B— 1 ) 


where  consideration  is  limited  here  to  the  purely  quadratic  case; 
see  bottom  of  page  8.  Combining  these  results  leads  to  FO  PD 


3//(  X) 
3X(  m) 


E{g'  C(m)  g  exp[g'  D(X)  g ] }  for  m=l:M  , 


( B-2  ) 


where  symmetric  matrix  D(X)  is  given  in  equation  (19)  as 


M 

D( X)  =  YU  Mm)  c(m)  •  (B-3) 

m=l 


Perform  the  same  eigen-decomposition  on  matrix  D(X)  as  in 
equation  (48);  namely,  D(X)  =  V(X)  E(X)  V(X)'  =  V  E(X)  V' ,  and 
define  the  linearly  transformed  Kxl  zero-mean  Gaussian  RV 

y  -  vf  g  with  E{y}  =  0  ,  E{y  y' }  =  E{V'  g  g'  V}  =  I  .  (B-4) 


Then,  using  g  = 


3 /#(  X)  = 
3X(  m) 


E{yr 


=  E  {y ' 


V  y,  the  FO  PD  in  equation  (B-2)  becomes 

V'  C(m)  V  y  exp(y'  V'  D(X)  V  y}  = 

F(X,m)  y  exp ( y '  E(X)  y) }  for  m=l:M  , 


( B-5 ) 


where  E(X)  =  diag{ek(X)},  and  symmetric  KxK  matrix 


B— 1 


F(X,m)  ■  V'  C(m)  V  =  V(X)'  C(m)  V(X)  for  m=l:M  . 


( B-6  ) 


Define  the  elements  of  this  KxK  F  matrix  as 

F(X,m)  =  [ f ( X,m;k,k ) ]  for  k,k=l:K  ;  m=l:M  ,  (B-7) 

and  let  the  components  of  RV  y  in  equation  (B-4)  be  denoted  as 
y  =  [y(l)  •••  y(K)]f  .  (b-8) 


Then,  the  FO  PD  in  equation  (B-5)  becomes,  for  m=l:M, 


H  -  e(  YZ  f  (  X,m;k ,k )  y(k)  y(k)  exp 
Ik ,  k=l 


r  K 


ZZ  eD(X)  y(p} 

lp=l  P 


TZ  f ( X,m; k , k )  E^y(k)  y(k)  exp 

k ,  k=l  l 


K 


1  I  e  (X)  y(p)' 
Ip-i  p 


( B-9  ) 


If  k  7*  k,  the  expectation  in  equation  (B-9)  is  zero.  Therefore, 


K  f  i\ 

=  ZZ  f  (  X,  m;k  ,  k  )  E<Jy(k)2  exp  ZZ  eD(x)  Y^P) 


k=l 


lp=l 


.  ( B-10 ) 


The  k-th  average  in  equation  (B-10)  is,  with  the  help  of  the 
statistics  of  Gaussian  RV  y  in  equation  (B-4), 

K  K 

(1  -  2  e,  )"3/2  T7  (1  -  2  e  )~h  =  (1  -  2  e,)-1  TT  d  -  2  e  )~H  = 
K  p=l  P  K  p=l  p 

p^k 

-  (1  -  2  ek)-1  /j(\)  for  k-l:K  .  (B-ll) 

This  last  result  for  joint  MGF  /j(\)  follows  from  equations  (23) 
and  (50)  in  the  purely  quadratic  case.  The  use  of  equation 


B-2 


( B— 11 )  in  equation  (B-10)  yields 

K 

ftfi}  "  0(X)  C  f ( X,m; k , k )  [1-2  e^X)]-1  for  m=l:M  .  (B-12) 


But,  since  X(X)  =  log  //(X), 

jXjX)  _  f ( X, m;k , k ) 

3X(m)  1-2  ek(X) 


there  immediately  follows 

for  m=l:M  . 


( B-13 ) 


The  elements  {f(X,m,k,k)J  of  matrix  F(X,m)  are  given  in  equations 
(B-6)  and  (B-7).  If  the  KxK  eigenvector  matrix  V(X)  in  equations 
(B-4)  through  (B-6)  is  expressed  in  terms  of  its  Kxl  column 
vectors  {Vk(X)},  k=l:K,  according  to  V(X)  =  (V^(X)  VR(X)], 

then  equations  (B-6)  and  (B-7)  yield 

f ( X,m, k , k )  =  Vk(X)r  C(m)  Vk(X)  for  m-lsM  ,  k=l:K  ,  (B-14) 

which  avoids  the  calculation  of  the  entire  KxK  F(X,m)  matrix  for 
each  m. 


The  result  in  equation  (B-13)  can  be  manipulated  into  a 

!' 

familiar  form: 

=  tr  { [  I  -  2  E  (  X )  ] -1  F(  X,m) }  =  tr{[I  -  2  E(X)]_1  V'  C(m)  V} 

=  tr { V  [1-2  E( X) ]-1  V'  C(m) }  =  tr {Q( X)-1  C(m)}  ,  (B-15) 

where  equations  (B-6),  (B-7),  and  (49)  were  used.  This  latter 
result  in  equation  (B-15)  agrees  with  equation  (52)  because  the 
second  part  of  the  joint  CGF,  X2<X),  is  zero  in  this  purely 
quadratic  case. 
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SECOND-ORDER  PARTIAL  DERIVATIVES  OF  CGF 

This  presentation  will  be  somewhat  abbreviated,  since  the 
details  are  similar  to  those  above.  Using  shorthand  notation 


Vx> 


9£jX) 

3X(m)  ' 


VX> 


3X(X)  //m(X) 

3X(m)  fj(\) 


( B-16 ) 


There  follows,  for  the  second-order  PDs  of  the  joint  CGF, 


//  (X) 

mm 


xmm(X)  ~  fj(\)  Xm(X)  Xm(X)  * 


( B— 17  ) 


From  equation  (46),  the  pertinent  quantity  is 


/t/mm(X)  =  exP  HH  X<P)  Z<P) 

—  I  lp=l 


kk  11 


=  E  S  f(n;k,k)  f  ( m;l ,  1 )  E<jyk  yk  yx  y1  exp|^j  ep  yp  (  '  (B-18) 


Only  two  cases  yield  nonzero  averages  in  equation  (B-18).  In 
the  first  case,  k  =  k  =  1  =  1^,  the  statistical  average  becomes 

E{y£  exp(g  6p  yp)}  '  E(yk  expK  yk])  TJ  E(exp[eP  yp])  - 

p^k 

K 

*  - 2 - E~ry  TT  [1  -  2  e.)~h  =  - 3 . £l±L—  for  k-l:K.  (B-19) 

[1-2  e.  ]3/Z  p-1  '  K  [1-2  e.p 

P^k  k 

For  convergence  of  these  integrals,  it  is  necessary  that  all  the 
eigenvalues  ek  =  ek(X)  <  H  for  k=l:K.  This  case  contributes  the 


B-4 


following  term  to 


K  f  ( m ;  k  ,  k )  f  (  m ;  k  ,  k  ) 

3  /j(\)  Y~[  - 5 -  for  •  (B-20) 

k=l  [1  -  2  ekr 


The  second  case  consists  of  three  subcases: 


k  =  k  /  1  =  1  , 
k  =  1  /  k  =  1  , 
k  =  1  /  k  =  1  . 


( B-21 ) 


The  pertinent  average  for  the  first  subcase  is 


2  2 

Elyk  yl  exp 


K 

,C  e 
lp=l 


p  yp 


fl(X) 


[1-2  ek][l  -  2  e1] 


( B-22 ) 


This  first  subcase  contributes  the  term 


K  f(m;k,k)  f ( m ; 1 , 1 ) 
"<X)  5  11  -  2  MU  -  2 


( B-23 ) 


=  fj{\) 


K 

r — i  nm 

[hi  1  -  2  eu 


f  ( m ;  k  ,  k  ) 


'  K  f  ( m ;  k ,  k ) 1 

1  —  2  e. 
U=1  k; 


K  f  ( m ;  k  ,  k )  f  ( m ;  k  ,  k  ) 

YZ  ■ 

k=l 


[1  -  2  ek] 


The  other  two  subcases  each  contribute  the  following  term  to 


/J  mm(X) 
mm 


K  f(m;k,k)  f(m;k,k) 
"<X)  S  [1  -  2  ekHl  -  2  ekl 


(B-24) 


(  K  f(m;k,k)  f(m;k,k) 

'  "tX)  \k£l  [1  -  2  vTTT  -  2  ek) 


K  f(m;k,k)  f(m;k,k)' 

C  - ~~~2 - 

k=l  [1-2  ekr 
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When  all  the  terms  are  combined  according  to  equation  (B-17), 
a  number  of  cancellations  occur,  yielding,  for  m,m=l:M, 

K  f(X,m;k,k)  f(X,m;k,k) 

Xmm(X)  =  2  [1-2  ek(X)][l  -  2  ek(X)]  #  (B-25) 

where  the  X  dependence  has  been  reintroduced,  and  equation  (B-13) 
was  used. 

Equation  (B-25)  may  be  manipulated  as  follows: 

X^(X)  =  2  tr  { ( I  -  2  E ) -1  F(X,m)  (1-2  E)"1  F(X,m))  = 

=  2  tr { ( I-  2  E)-1  V'  C ( m)  V  (I  -  2  E)-1  V'  C(m)  V}  ,  (B-26) 

upon  use  of  equations  (B-6)  and  (B-7).  Then,  upon  movement  of 
the  trailing  matrix  V  to  the  front  of  the  trace,  there  follows 

X  (X)  =  2  tr{Q( X)-1  C ( m)  Q(X)-1  C(m)}  for  m,m=l:M  ,  (B-27) 

mm  ~  ■“ 

where  equation  (49)  was  used.  This  result  is  identical  to 
equation  (56)  in  this  purely  quadratic  case. 


RELATED  EIGENVALUE  PROPERTIES 

Suppose  matrix  D  is  KxK  and  symmetric.  Its  eigen- 
decomposition  is  D  =  V  E  V'  or  D  ^  Vk  for  k=l:K.  Then, 

V£  D  =  ek  V£  and  V1  =  &kl.  Equivalently,  eR  ■  V£  D  . 

Now,  suppose  that  matrix  D  is  a  function  of  scalar  x.  Then, 
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Relation  (B-29)  is  true  for  any  matrix  D(x).  Now,  let 
X  =  [ X( 1 )  • • •  X(M) ] '  and 


M 

D  =  D( X)  =  YZ1  Mm)  C(m)  /  KxK  C(m)  constant  ,  (B-30) 

m=l 


and  interpret  x  as  X(m).  Then,  D(X)  Vk(X)  =  ek(X)  Vk(X)  and 


3ek(X) 
3X(  m) 


VX>' 


3D(  X) 
3X(m) 


Vk(X)  =  Vfc(X)'  C(m)  Vk(X) 


( B-31 ) 


for  k=l:K,  m=l:M.  This  is  a  useful  relation  for  the  PD  of  an 
eigenvalue.  The  quantity  on  the  right-hand  side  of  equation 
(B-31)  is  identical  to  the  quantity  in  equation  (B-14). 


Upon  multiplication  of  equation  (B-31)  by  X(m),  summation 
over  m,  and  use  of  equation  (B-30),  there  follows 

M  3ek(X) 

HJ  Mm)  m)  =  Vk(X)'  D(  X)  Vk(X)  =  ek(X)  for  k=l:K  .  (B-32) 
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This  result  gives  an  expansion  of  an  eigenvalue  in  terms  of  its 
PDs  weighted  by  the  {X(m)}. 


A  second  PD  on  equation  (B-31)  yields 


9  ek(X) 


9Vk(X) ' 


9Vk(X) 


axim)  3X(n)  -  8X(n)  C(in|  VX)  +  VX>'  C"”>  8X? ~  '  IB-33) 


Multiplication  by  X(m)  and  summation  over  m  now  leads  to 


M  92ek(X)  9V,  (X)'  9V,(X) 

C  Mm)  yxoa)  8X(n)  =  IXuT  D(X)  VX>  +  VX>'  D(X>  alfe" 


9Vk(X) ' 


9Vk(X) 


=  -axur  Vx>  ek<x>  +  ek(X)  Vx>'  IxfnT 


=  ek(X)  3X(n) [Vk(X) '  Vk(X)]  =  0  for  k=1:K  »■  n=l :M  .  (B-34) 


Thus,  the  sum  of  second-order  PDs  of  the  eigenvalues,  weighted  by 
the  { X(m) } ,  is  always  identically  zero.  Relations  (B-31), 

(B-32),  and  (B-34)  have  been  checked  numerically. 
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APPENDIX  C  -  MATRIX  PROPERTIES 


TRACE  PROPERTY 

Let  A  be  a  KxK  matrix,  not  necessarily  symmetric.  The  trace 
of  matrix  A,  in  terms  of  its  elements  {A(k,l)},  is 

K 

tr ( A)  =  XI]  A(k,k)  .  ( C-l ) 

k=l 

Then,  the  trace  of  a  product  of  two  KxK  matrices  follows  as 

K 

tr ( A  B)  =  XU  A(k,l)  B( 1 , k )  =  tr(B  A)  .  ( C-2 ) 

k ,  1=1 

It  immediately  follows  that  the  trace  of  a  product  of  several 
KxK  matrices,  such  as  C  D  E  F  G,  can  be  rearranged,  for  example, 
as 

tr ( C  D  E  F  G)  =  tr ( E  F  G  C  D)  ,  (C-3) 

simply  by  identifying  A  as  C  D  and  identifying  B  as  E  F  G.  In 
fact,  any  cyclic  rearrangement  of  the  matrices  is  allowed  without 
changing  the  value  of  the  trace.  However,  if  the  cyclic  pattern 
is  changed  (for  example,  by  switching  the  locations  of  matrices  F 
and  G),  the  trace  is  modified. 
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EIGENVALUE  PROPERTY 


Let  matrix  A  be  KxK  and  have  eigen-decomposition 

A  =  V  E  V-1  ,  E  =  diag{ek}  ,  k=l:K  .  (C-4) 

Then,  A2  =  A  A  =  V  E  V-1  V  E  V-1  =  V  E2  V-1,  which  can  be 
immediately  generalized  to 

AP  =  V  EP  V-^  for  integer  p  .  (C-5) 

There  follows 

K 

tr(AP)  =  tr ( V  EP  V-1)  «  tr(EP  V-1  V)  =  tr(EP)  =  yf]  eP  ( C-6 ) 

k=l  K 

Thus,  the  trace  of  the  p-th  power  of  A  is  equal  to  the  sum  of  the 
p-th  powers  of  all  the  eigenvalues  of  matrix  A. 

USEFUL  MATRIX  PROPERTIES 

Let  A  be  a  KxK  matrix  as  in  equation  (C-4).  For  scalar  a, 
the  identity 

(1  -  a)-1  =  a  =  1  +  a  +  a2  +  a3  +  • • •  if  j a |  <  1  .  (C-7) 

By  using  equation  (C-4)  and  I-A-V(I-E)  V-1,  this 
generalizes  to  the  matrix  relation 

(I  -  A)-1  =  V  (I  -  E ) -1  V-1  =  v  diag{ ( 1  -  eR ) -1 >  V-1  = 

=  V  diag{ 1  +  ek  +  ©k  +  • • • }  V-1  =V(I+E+E2+***}  V-1  = 

=  I  +  A  +  A2  +  A3  +  • • •  if  | ek |  <  1  for  k=l:K  .  (C-8) 
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That  is, 


C  An  = 
n=0 

where  eig(A) 


( I  -  A)  1  if  |eig(A) |  <  1 
denotes  all  the  eigenvalues 


of  matrix  A. 


( C-9 ) 


A  function  f(A)  of  matrix  A  is  defined  according  to 

f ( A)  s  V  f(E)  V-1  where  A  =  V  E  V-1  .  (C-10) 

2  3 

Therefore,  the  relation  -log(l  -  a)  =  a  +  a  /2  +  a  /3  +  •••  for 
scalar  a  generalizes  to 

-log ( I  -A)  =  -V  log( I  -  E)  V-1  =  -V  diag{log(l  -  ek)}  V-1  = 


TT  ..  (  12  1 
=  -V  ciiag[-ek  -  jek  -  ye 


=  A  +  yA3  +  yA3  +  •••  if  |ek|  <  1  for  k=l:K  . 


(C-ll) 


That  is, 


5  !  -An  =  -  log( I  -  A)  if  |eig(A)|  <  1  .  (C-12) 

n=l 


Now,  let  A( x )  be  a  KxK  matrix  which  is  a  function  of  x. 
Represent  the  derivative  with  respect  to  x  by  the  symbol  A(x). 
Then,  by  use  of  the  chain  rule,  the  derivative  of  the  p-th  power 
of  A(x)  contains  p  terms: 

af  A(X)P  -  af(A(x)  •••  A(x|)  =  A  aP‘1  +  a  a  aP‘2  +  •••  +  A'”1  A  • 

( C-13 ) 
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Then,  by  using  the  trace  properties  in  equations  ( C-2 )  and  (C-3), 
there  follows 

tr^-A(x)PJ  =  p  tr^A(x)p_1  A(  x  )  j  .  (  C-14  ) 

If  matrix  B(x)  is  the  inverse  of  matrix  A(x),  the  derivative 
of  matrix  B(x)  may  be  found  as  follows: 

B ( x )  -  A(x)"1  ,  B ( x )  A( x )  =  I  ,  B ( x )  A( x )  +  B(x)  A(x)  =  0  , 
B(x)  =  -  B(x)  A(x)  A(x)  ^  =  -  B ( x )  A(x)  B(x)  .  (C-15) 

That  is,  the  derivative  of  the  inverse  of  a  matrix  A(x)  is  given 
by  the  negative  derivative  of  the  matrix  A(x),  which  is  then  pre- 
and  post-multiplied  by  the  inverse  matrix  B(x). 

The  final  needed  matrix  property  involves  the  derivative  of 
equation  (C-ll).  There  follows 

gf[-log(I  -  A)]  -  A  +  i(A  A  +  A  a)  +  A2  +  A  A  A  +  A2  a)  +  ... 

( C-16 ) 

from  which  is  obtained,  by  using  equation  ( C— 8 ) , 

tr^-iogd  -  A)])  -  tr(A  ♦  A  A  +  A2  A  +  •••)  -  tr((I  -  A)-1  a). 

( C— 17 ) 

Finally,  interchanging  the  trace  and  derivative, 

^[trf-logll  -  A(x)]]]  =  tr([I  -  A(x)]-1  A(x)}  .  (C-18) 

This  holds  for  all  A(x),  except  if  one  or  more  eig(A(x))  =  1. 
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APPENDIX  D  -  MATLAB  PROGRAM  quadlinspa 


clear  al 1  %  SPA  tc 

i  joint  PDF  of  M  quadratic  and  linear 

]yf=4;  %  Number  of  quadratic  and  linear  forms. 

K=64;  %  Number  of  Gaussian  random  variables . 

tol=le-7;  %  Tolerance  in  saddlepoint  search. 

kkrrax=100;  %  Maximum  number  of  search  trials. 

f=.499;  %  Proximity  to  boundary  at  .5 

randn( 'state' ,0) 

%  USIEUr  INPORMAITCN 

A=randn(K,K)  ; 

%  Positive-definite  covariance 

R=A*A' ; 

%  matrix,  R,  of  K  Gaussian  FVs. 

r=randn(K,  1)  ; 

%  Mean  vector,  r,  of  K  Gaussian  RVs. 

P=zeros  (K,K,M) ; 

for  rrt=l:M 

A=randn(K,K)  ; 

P ( : , :  ,m)  =  (A+A1 )  * . 5; 

%  symmetric  quadratic  terms,  P 

erd 

p=randn(K,M)  ; 

%  Linear  terms,  p 

q=randn(M, 1)  ; 

%  Constant  terms,  q 

z=zeros  (M,  1)  ; 

%  SPECIFY  FIELD  POINT  z 

S=chol (R)  ; 

%  KxK 

g=randn(K,l)  ; 

%  Kxl,  N(0,1) 

w=r+S '  *g; 

%  Kxl,  N(r,R) 

for  m=l:M 

z(m)  =  (P(: , :  ,m)*w+-p( 

:  ,m) )  *w+q(m)  ; 

erd 

%  Mxl,  field  point  z 

C=zeros  (K,K,M)  ; 

%  ERE-CCMFJIATICN  OF  MAIRICES 

v-zeros  (K,M)  ; 

c=zeros  (M,  1)  ; 

S=chol  (R)  ; 

%  KxK 

for  rrt=l:M 

A=P ( : ,  :  ,m) ; 

%  KxK 

C  ( : ,  :  ,m)  =S*A*S ' ;  • 

%  K«C 

a=A*r; 

%  Kxl 

v(:,m)=S*a; 

%  Kxl 

c(m)=r'*a; 

%  1x1 

aid 

v=S*p+2*v; 

%  K>M 

c=q+p'*r+c; 

%  F&d 

tic 

%  SEARCH  FOR  SADDLEPOINT 

L=zeros  (M,  1) ; 

B=zeros(K,K,M)  ; 

G=zeros  (M,  1)  ; 

H=zeros  (M,M)  ; 
vt=v'  ; 

%  MxK 

kk=0; 

K2=K*K; 

znorm=sqrt  (z '  *z)  ; 
err=z-G; 

%  MdL 

while  (sqrt  (err '  *err)  /znorm)  >tol 
t=v*L; 

%  Kxl 

P=reshape(C,K2,M)  ; 

DL=resbape  (P*L,  K,  K)  ; 

%  KxK 

e=eig(DL)  ; 

%  KxL 

erraBx(e)  ; 
if  src>=.5 

L=L*(f/en)  ; 

%  Mxl 

DD=DL*  (f/ena)  ; 

%  KxK 

eigmax=[an  kk] 
e±l 

Q=eye(K)-2*DL; 

%  KxK 

qt=Q\t; 

%  Kxl 

for  ro=l:M 

B1=C ( : , : ,m)  ; 

%  KxK 

A=Q\B1; 

%  KxK 

B  ( : , :  ,rti)  =A; 

G  (m)  =trace  (A)  +gt '  *Bl*qt ; 
end 

G=G+vt*at+c;  %  Mxl  Gradient  vectoi 

for  ml=l:M 

B1=B ( : , :  ,ml) ; 

%  KxK 

tb=2*t  ‘  *Bl+v( :  ,ml) 

%  IxK 

for  m2=ml:M 

E2=B( : , :  ,m2) ; 

%  KxK 

ts=Bl ( : ) ' ^reshape (B2 

: ' ,K2,1) . . . 

+  (tb*B2+v(:  ,m2)  '*B1) 

*qt; 

H(ml,m2)=ts; 

H(rn2,ml)=ts; 

erl 

exl 

H=H*2+vt/Q*v;  %  KfcM  Hessian  matrix 

err^z-G; 

%  McL 

dI>H\err; 

%  Md 

fr=.3;  %  fraction:  [0  1) 
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%  Md 


f  f=l-fd (kk+1) ; 

L=b+dL*ff; 
kk=kk+l; 

if  kk>kkmax,  break,  end 
end  %  while 

disp  ( [ '  kk  =  ’int2str(kk)  ] ) 

L=L+dL*(l-ff) ;  %  saddl^oint  %  MxfL 

u=c 1  *L;  %  bd 

t=V*L;  %  Kxl 

P=resbape(C,K2,M)  ; 

Dl^reshape  (P*L,  K,  K) ;  %  KxK 

e=eig(DL)  ; 
if  (max  (e)  >f) 

disp ( [ ' eigmax  is  greater  than  f  =  ’num2str(f)  ] ). 
keyboard 
exl 

Q=eye(K)-2*DL; 
qt=Q\t; 

br=.5*(t'*qt)+u; 
rrgfO=l/sgrt  (prod(l-2*e) )  ; 
rrigf=ingf  0*exp  (hr)  ; 
cgf =log  (mgf  0 )  +br  ; 

for  m=l:M 

Bl=C(:,:,m); 

A=Q\B1; 

B  ( : , :  ,m)=A; 

G  (m)  =trace  (A)  +qt '  *Bl*qt; 
end  , 

G=G+vt*gt+c;  %  MxL  Gradient  vector. 

err=z-G;  %  Error  in  gradient  of  CGF. 

reg=sqrt  (err '  *err)  /znorm; 
disp  ( [ '  rel_err_grad  =  '  num2str  (reg)  ] ) 

tl=tOC; 

disp( [ 'tl(sec)  =  'num2str (tl) ] ) 
tic 

BB=zenos  (K,K,M,M)  ; 
for  ml=l:M 

B1=B( : , : ,ml) ;  %  KxK 

tb=2*t '  *Bl+v( :  ,ml) ' ;  %  biC 

for  m2=l:M 


%  KxK 
%  KxL 
%  bd 


%  KxK 
%  KxK 


B2=B(: , : ,m2) ; 

%  KxK 

A=B1*B2 ; 

%  KxK 

BB(; ,  :,ml,m2)=A; 

if  (ml<=m2 ) 

ts=trace  (A) . . . 

%  1x1 

+  (tb*B2+v  ( :  ,m2 )  '  *B1)  *qt; 

H(ml,m2)=ts; 

H(iri2,ml)=ts; 

aid 

aid 

aid 

qv=Q\v;  %  KxM  • 

H=H*2+vt*qv;  %  MxM  Hessian  matrix 

den=sqrt  ( (2*pi)  ^M*det  (H) )  ; 
pdfO=mgf*e^(-z'*L)/den;  %  SPAO 

T=zeros  (M,M,M)  ; 
for  ml=l:M 
for  m2=ml:M 

A=reshape  (BB  ( : :  ,ml,m2 )  ' , 1, K2 ) ; 
for  m3^n2  :M 

T(ml,ni2,mB)=A*reshapefB(: , :  ,m3)  ,K2,1)  ; 

end,  end,  end 

for  ml=l:M 

for  m2=l:M 

for  m3=l:M 

s=sort  ( [ml  m2  m3]); 

T(ml,m2,m3)=T(s (1)  ,s(2)  ,s(3) )  ,- 
end,  end,  end 

1>T*8;  %  MxMxM;  Third-Order  Partial  Derivatives 

Tl=zeros  (M,M,M)  ; 

T2=zeros  (M,M,M)  ; 

T3=zeros  (M,M,M)  ,- 
for  m=l:M 

A2=qv'*C(:,  :,m)*qv;  %  M&L 

Tl(m,  :)=A2; 

T2  ( :  ,m,  :  )=A2; 

T3  ( : , :  ,m)=A2; 
ed 

Th=  (T1+T2+T3 )  *2 ;  %  iyfcMxM;  TO  EDs 


D-4 


for  ml=l:M 
for  m2=l:M 

Al=vt*(  (BB(: ,  :  ,ml,rn2)+BB( : , :  ,m2,ml)  )*qt)  ; 
Tl(:,ml,m2)=Al;  %  Md. 

12  (m2, :  ,ml)=Al; 

13  (ml, m2,  :  )=A1; 
end,  end 

1b=  (H+12+13 )  *4 ;  %  10  EDs 

for  ml=l:M 
for  m2=l:M 

Bl=t'*EB( : , :  ,ml,m2) ;  %  IxK 

for  m3=l:M 

a=Bl*B(: , :  ,m3)  *gt;  %  lxl 

11  (ml,  m2,  m3)  =a; 

12  (m3,ml,m2)=a; 

13  (m2, m3, ml)  =a; 
end,  end,  end 

Tt=  (11+12+13)  *8;  %  10  EDs 

T=T+1d+1b+Tc;  %  MxMxM;  Ihird-Qrder  Partial  Derivatives 

F=zeros(M,M,M,M)  ; 
for  ml=l:M 
for  m2=ml:M 

A=reshape (BB ( : , :  ,ml,m2)  '  ,1,K2)  ; 
for  m3=m2:M 

Bl=reshape (BB ( : , :  ,ml,m3)  '  ,1,K2)  ; 
for  m4=m3:M 

F(ml,m2,m3,m4)=A*reshape(BB( : , :  ,m3,m4) . . . 

+BB( : , :  ,m4,m3)  ,K2,1)  . . . 

+Bl*reshape (BB ( : , :  ,m2,m4)  ,K2,1)  ; 

end,  end,  end,  end 

for  ml=l:M 

for  m2=l:M 

for  m3=l:M 

for  m4=l:M 

s=sort  ( [ml  m2  m3  mi] )  ; 

F(ml,m2,m3,m4)=F(s (1)  ,s(2)  ,s(3)  ,s(4) )  ; 
end,  end,  end,  end 

F=F*16;  %  Fouirth-Order  Partial  Derivatives 

H=zeros  (M,M,M,M)  ; 
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12=H;  T3=T1;  T4=H;  15=H;  T6=H;  T7=H; 

T8=n;  T9=H;  T10=n;  T11=!E1;  H2=H; 

for  ml=l:M 
for  m2=l:M 

A2=vt* (BB(: ,  :,ml,m2)+EB(:,  :,m2,ml))*qv-; 

T1  (ml, m2 , : , : )  =A2;  %  E&M 

T2  (ml, :  ,m2 , : )  =A2  ; 

T3(ml, : , :  ,m2)=A2; 

T4(:  ,ml,m2, :)=A2; 

15 (:  ,ml, :  ,m2)=A2; 

T6(: ,  :,ml,m2)=A2; 
end,  end 

1^=(11+T2+13+T4+T5+T6)*4;  PO  EDs 

for  ml=l:M 
for  m2=l:M 

Bl=vt*(BB(:, :  ,ml,m2)+BB(: , :  ,m2,ml) ) ;  %  MxK 
for  m3=l:M 

A1=B1*  (B( : , :  ,m3)  *qt) ;  %  Mxl 

H(:  ,ml,m2,m3)=Al; 

12  (m3, :  ,ml,m2)=Al; 

13(m2,m3,:,ml)=Al; 

T4 (ml, m2, m3,  :)=A1; 

T5(:,ml,m3,m2)=Al; 

T6  (m2 , :  ,ml ,m3 )  =A1  ; 

T7  (m3, m2, :  ,ml)=Al; 

T8 (ml, m3, m2,  :)=A1; 

T9  (m3, m2, ml,  :)=A1; 

HO  (m2, ml, :  ,m3)=Al; 

HI  (ml, :  ,m3,m2)=Al; 

H2  ( :  ,m3,m2,ml)=Al; 

end,  end,  end 

Tb=  (H+12+13+T4+T5+T6+ . . . 

T7+T8+T9+H0+H1+H2)  *8;  %  MxI>M£<M;  PO  PDs 

for  ml=l:M 
for  m2=l:M 

Bl=t '*  (BB(: , :  ,ml,m2)+BB( : , :  ,m2,ml) ) ;  %  lxK 
for  m3=l:M 
for  m4=l:M 

a=Bl*BB  ( : , :  ,m3  ,m4)  *qt;  %  1x1 

H  (ml, m2  ,m3  ,m4)  =a,- 
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T2  (m4,ml,m2,m3)=a; 

T3  (m3,m4,ml,m2)=a; 

T4(m2,it8,nv4,inl)=a; 

T5(ml,m2,m4,m3)=a; 

T6(ir6,inl,rn2,rn4)=a; 

T7  (m4r,m3  ,rrLL,xn2 )  =a.; 

T8  (rtt2  ,m4rm3  ,ml)  =a; 

T9(ml,m3,in2,m4)=a; 

H0(m4,ml,m3,iTi2)=a; 

HI  (m2  )  =a,- 

T12  (m3  ,m2  ,m4,ml)  =a; 
end,  end,  end,  end 
Tc=  (H+T2+T3+T4+T5+T6+. . . 

H+T8+T9+H0+T11+T12)  *8;  %  FO  EDs 

F=F+TS.+'Ii>f Tt ;  %  MMM;  Fourth-Order  Partial  Derivatives 
%  CALCULATE  ODEPBCITCN  TEEMS 

A2=zeros  (M,M)  ,- 
M2=M*M; 

Hi=inv(H) ;  %  ^ 

Hr=Hi  ( : )  ' ;  %  1*^2 

for  ml=l:M 
for  m2=l:M 

A2  (ml,m2)=Hr*reshspe(F(: , :  ,ml,m2)  ,M2,1)  ; 
end,  end 
c4=Hr*A2  ( : )  /8; 

Al=zeros  (M,  1)  ; 
for  itt=l:M 

A1  (m)  =Hr*reshape  (T  ( : , :  ,m)  ,M2 , 1)  ; 
erl 

c3a=-Al'*Hi*Al/8; 

A3=zeros  (M,M,M)  ; 
for  ml=l:M 
B2=Hi( :  ,ml) '  ; 
for  nn2=l:M 
B3=Hi(:  ,m2)  ; 
for  m3=l:M 

A3  (ml,m2,m3)=B2*T( : , :  ,m3)*B3; 
end,  end,  end 
B2=zeros  (M,M)  ,- 


%  Mxl 

%  mi 
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for  ml=l:M 

B3=reshape(T( : , :  ,ml)  ,1,M2) ; 
for  m2=l:M 

B2  (ml ,m2 )  =B3*reshape (A3  ( : , :  ,m2)  ,M2,1) ; 
end,  end 

c3b=-Hr*B2  ( : )  /12  ; 

ct=c4+c3a+c3b;  %  FTRST-ORDER  CORBBCZECN  TERM 

disp  ( '  c4  c3a  c3b  ct  = 1 ) 

disp([c4  c3a  c3b  ct]) 

pdfl=pdf 0* (1+ct) ; 

pdfe=pdfO*exp(ct) ; 

pdfb^pdfO*  (l+ct/2)  /  (l-ct/2) ; 

t2=toc; 

disp  (['t2 (sec)  =  'num2str(t2) ] ) 
disp(  ['pdfO  =  'num2str(pdf0)  ] ) 
disp(['pdfl  =  ‘num2str(pdfl)  ] ) 
disp  ([ 'pdf  e  =  'num2str(pdfe)  ] ) 
disp(['pdfb  =  ’nurn2str(pdfb)  ] ) 
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APPENDIX  E  -  MATLAB  PROGRAM  quadspa 


clear  all 
M=4; 

K=7; 

tol=le-7; 
kkrrax=100  ; 
f=.499; 


%  SPA  to  joint  PDF  of  M  quadratic  forms. 
%  Number  of  quadratic  forms. 

%  Number  of  Gaussian  random  variables . 

%  Tolerance  in  saddlepoint  search. 

%  Maximum  number  of  search  trials . 

%  Proximity  to  boundary  at  .5 


randn(  'state'  ,0) 
A=randn(K,K)  ; 

R=A*A' ; 

P=zeros  (K,K,M)  ; 
for  rrt=l:M 
A=randn(K,K) ; 

P  ( : ,  :  ,m)  =  (A+A ' )  * .  5 ; 
aud 


%  INFUT  INFORMATION  R  AND  P 
%  Positive-definite  covariance 
%  matrix,  R,  of  K  Gaussian  RVs. 


%  Svrrmetric  quadratic  terms,  P 


z=zeros  (M,  1)  ; 

S=chol (R)  ; 
g=randn(K,  1)  ; 

S  ^ 
for  m=l:M 

Z  (m)  =w'  *P  ( : , :  ,m)  *W; 
erd 


%  SPECIFY  FIELD  PDINT  z 
%  K^C 

%  Kxl,  N(0, 1) 

%  Kxl,  N(0,R) 


%  Mxl,  field  point  z 


C=zeros(K,K,M) ;  %  PRE-GCMPUTAITCN  OF  MATRIX  C 

S=cbol  (R) ;  %  KxK 

for  m=l:M 

C( : ,  :  ,m)=S*P( : ,  :  ,m)  *S' ;  %  RxK 

ead 


tic  %  SEARCH  FOR  SADDLEPOINT 

L=zeros  (M,  1) ; 

B=zeros  (K,K,M)  ; 

G=zeros  (M,  1)  ; 

H=zeros  (M,M) ; 

Wc=0; 

K2=K*K; 

znorm=sqrt  (z '  *z)  ; 

err=z-G;  %  Mxl 

viiile  (sqrt  (err '  *err)  /znorm)  >tol 
P=reshape  (C, K2  ,M)  ; 

.  DL=reshape  (P*L ,  K,  K) ;  %  KxK 

e=eig  (DL) ;  %  Kxl 


enraax(e)  ; 
if  en>=.5 


L=L*(f/en)  ; 

%  Mxl 

DL=DL*  (f/en) ; 

%  KxK 

eigcrax=  [en.  kk] 

end 

Q=eye(K)-2*DL; 

%  KxK 

for  m=l:M 

A=Q\C  ( : , :  ,m) ; 

%  KxK 

B  ( : , :  ,m)  =A; 

G(m)=trace(A) ; 

%  Mxl  Gradient 

end 


for  ml=l:M 

A=reshape (B ( : , :  ,ml)  '  ,1,K2)  ; 
for  m2=mL:M 

ts=A*reshape (B( : , :  ,m2)  ,¥2,1)  ; 
H(ml,m2)=ts; 


H(m2,ml)=ts; 

end 


end 
H=H*2  ; 
err=z-G; 
dL=H\err; 
fr=.6; 

ff=l-fr~  (kk+1) ; 
L=L+dL*ff ; 


%  M>$[  Hessian  matrix 
%  ly&d 
%  b&d 

%  fraction:  [0  1) 

%  MxL 


kk=kk+l; 

if  kk>kknax,  break,  end 
end  %  while 


disp(['kk  =  '  int2str  (kk)  ] ) 


L=L+dL*  (1-ff) ;  %  saddlepoint  %  Mxl 

P=resbape(C,K2,M)  ; 

DL=resbape  ( P*L ,  K,  K) ;  %  KxK 

e=eig(DL)  ; 
if  (max(e)>f) 

disp ( [ ' eigrrex  is  greater  than  f  =  'num2str(f)  ] ) 
keyboard 
end 

Q=eye  (K)  -2*DL;  %  Kx K 

mgf=l/sgrt  (prod(l-2*e) )  ; 
cgf=log(mgf); 
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%  KxK 


for  m=l:M 

A=Q\C  ( : , :  ,m)  ; 

B  ( : , :  ,Itl)  =A; 

G(m)=trace(A) ;  %  Mxl  Gradient  vector 

end 

err=z-G;  %  Error  in  gradient  of  cgf 

reg=sgrt  ( err 1  *err )  /znorm; 

disp  ( [  ’  rel_err_grad  =  '  num2str  (reg)  ] ) 

tl=tOC; 

disp( [ 'tl(sec)  =  'num2str(tl)]) 
tic 

BB=zeros(K,K,M,M) ;  • 
for  ml=l:M 


B1=B( : ,  :  ,ml) ; 

%  KxK 

for  m2=l:M 

A=Bl*B(:,:,m2); 

%  KxK 

BB  ( : ,  :  ,ml,m2)=A; 

if  (ml<=m2) 

ts=trace(A) ; 

%  lxl 

H(ml,m2)=ts; 

H(m2,ml)=ts; 

end 

end 

end 

H=H*2;  %  Hessian  matrix 

den=sart  ( (2*pi)  ^M*det  (H) )  ; 

pafO=mgf *exp (-z '  *L)  /den;  %  SPAO 

T=zeros  (M,M,M)  ; 
for  ml=l:M 
for  m2=ml:M 

A=reshape (BB ( : , :  ,ml,ni2)  '  ,1,K2)  ; 
for  m3=m2  :M 

T(ml,m2 ,m3 )  =A*reshape (B ( : , :  ,m3)  ,K2,1)  ; 

end,  end,  end 

for  ml=l:M 

for  m2=l:M 

for  m3=l:M 

s=sort ( [ml  m2  m3]); 

T(ml,m2,m3)=T(s(l)  ,s(2)  ,s(3) )  ; 
end,  end,  end 


T=T*8;  %  MxM&dM;  Third-Order  Partial  Derivatives 

F=zeros  (M,M,M,M)  ; 
for  ml=l:M 
for  m2=ml:M 

A=reshape  (BB  ( : , :  ,ml,m2)  '  ,1,K2)  ; 
for  m3=5n2  :M 

Bl=reshape (BB ( : , :  ,ml,m3)  '  ,1,K2)  ; 
for  m4=m3  :M 

F(ml,m2,m3,m4)=A*reshape(BB(: ,  :,m3,m4) . . . 

+BB( : , :  ,m4,m3)  ,K2,1)  . . . 

+Bl*reshape (BB ( : , :  ,rti2 ,m4)  ,K2 , 1) ; 

end,  end,  end,  end 

for  ml=l:M 

for  m2=l:M 

for  m3=l:M 

for  m4=l:M 

s=sort  ( [ml  m2  m3  m4] )  ,- 
F(ml,m2,m3,m4)=F(s(l)  ,s(2)  ,s(3)  ,s(4) )  ,- 
end,  end,  end,  end 

F=F*16;%  MxMxffoM;  Fourth-Order  Partial  Derivatives 

%  CAIDXAIE  CORRECTION  TERMS 
A2=zeros  (M,M)  ; 

M2=M*M; 

Hi=irrv(H) ;  %  Mtf-1 

Hr=Hi  ( : )  ' ;  %  1*M2 

for  ml=l:M 
for  m2=l:M 

A2  (ml,m2)  =Hr*reshape(F( : , :  ,ml,m2)  ,M2,1); 
end,  end 
c4=Hr*A2  ( : )  /S; 


Al=zeros  (M,  1)  ,- 
for  m=l:M 

A1  (m)  =Hr*reshape (T(: , :  ,m)  ,M2,1)  ; 
exl 

c3a=-Al ' *Hi*Al/8 ; 

A3=zeros  (M,M,M)  ; 
for  ml=l:M 

B2=Hi(:,ml) %  Mxl 
for  m2=l:M 
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B3=Hi(:,m2);  %  M& 

for  m3=l:M 

A3  (ml, m2  ,m3 )  =B2*T( : , :  ,m3 )  *B3  ; 
end,  end,  end 
B2=zeros  (M,M)  ; 
for  ml=l:M 

B3=reshape(T( : , :  ,ml)  ,1,M2)  ; 
for  in2=l:M 

B2 (ml ,m2 )  =B3*reshape (A3 ( : ,  :,m2)  ,M2,1)  ; 
end,  end 

c3b=-Hr*B2(:)/12; 

ct=c4+c3a+c3b;  %  FTRST-OBDER  COHRECTICW  TERM 
disp ( ' c4  c3a  c3b  ct  =') 
disp(tc4  c3a  c3b  ct] ) 
pdfl=pdfO* (1+ct) ; 

pdfe=pdfO*ej^)(ct)  ; 
pdfb=pdfO* (l+ct/2) / (l-ct/2) ; 
t2=tOC; 

disp( [ 't2 (sec)  =  ’num2str(t2) ]) 
disp ( ['pdf 0  =  'num2str (pdfO)  ] ) 
disp ( ['pdf  1  =  'num2str(pdfl)  ] ) 
disp ([' pdf e  =  'num2str(pdfe)  ] ) 
disp(  [  'pdfb  =  'num2str(pdfb)  ] ) 
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